// 8-point/16-point windowed-sinc interpolation LUT generator

#include <stdint.h>
#include <stdbool.h>
#include <stdlib.h>
#include <math.h>
#include "ft2_windowed_sinc.h"
#include "../ft2_header.h" // PI
#include "../ft2_video.h" // showErrorMsgBox()

typedef struct
{
	double kaiserBeta, sincCutoff;
} sincKernel_t;

// globalized
float *fSinc[SINC_KERNELS], *fSinc8[SINC_KERNELS], *fSinc16[SINC_KERNELS];
uint64_t sincRatio1, sincRatio2;

static sincKernel_t sincKernelConfig[2][SINC_KERNELS] =
{
	/* Some notes on the Kaiser-Bessel beta parameter:
	** Lower beta = less treble cut off, more aliasing (narrower mainlobe, stronger sidelobe)
	** Higher beta = more treble cut off, less aliasing (wider mainlobe, weaker sidelobe)
	*/
	{ // -- settings for 8-point sinc --
		// beta, cutoff
		{  9.20, 1.000 }, // kernel #1 (beta < ~9.2 leads to audible aliasing here)
		{  8.50, 0.750 }, // kernel #2
		{  7.30, 0.425 }  // kernel #3
	},

	{ // -- settings for 16-point sinc --
		// beta, cutoff
		{  8.61, 1.000 }, // kernel #1 (beta 8.61 = Blackman-window approximation)
		{  8.50, 0.750 }, // kernel #2
		{  7.30, 0.425 }  // kernel #3
	}
};

// zeroth-order modified Bessel function of the first kind (series approximation)
static inline double besselI0(double z)
{
	double s = 1.0, ds = 1.0, d = 2.0;
	const double zz = z * z;

	do
	{
		ds *= zz / (d * d);
		s += ds;
		d += 2.0;
	}
	while (ds > s*(1E-12));

	return s;
}

static inline double sinc(double x, double cutoff)
{
	if (x == 0.0)
	{
		return cutoff;
	}
	else
	{
		x *= PI;
		return sin(cutoff * x) / x;
	}
}

// note: numPoints/numPhases must be 2^n!
static void makeSincKernel(float *out, int32_t numPoints, int32_t numPhases, double beta, double cutoff)
{
	const int32_t kernelLen = numPhases * numPoints;
	const int32_t pointBits = (int32_t)log2(numPoints);
	const int32_t pointMask = numPoints - 1;
	const int32_t centerPoint = (numPoints / 2) - 1;
	const double besselI0Beta = 1.0 / besselI0(beta);
	const double phaseMul = 1.0 / numPhases;
	const double xMul = 1.0 / (numPoints / 2);

	for (int32_t i = 0; i < kernelLen; i++)
	{
		const double x = ((i & pointMask) - centerPoint) - ((i >> pointBits) * phaseMul);

		// Kaiser-Bessel window
		const double n = x * xMul;
		const double window = besselI0(beta * sqrt(1.0 - n * n)) * besselI0Beta;

		out[i] = (float)(sinc(x, cutoff) * window);
	}
}

bool setupWindowedSincTables(void)
{
	sincKernel_t *k;
	for (int32_t i = 0; i < SINC_KERNELS; i++, k++)
	{
		fSinc8[i]  = (float *)malloc( 8 * SINC_PHASES * sizeof (float));
		fSinc16[i] = (float *)malloc(16 * SINC_PHASES * sizeof (float));

		if (fSinc8[i] == NULL || fSinc16[i] == NULL)
		{
			showErrorMsgBox("Not enough memory!");
			return false;
		}

		k = &sincKernelConfig[0][i];
		makeSincKernel(fSinc8[i], 8, SINC_PHASES, k->kaiserBeta, k->sincCutoff);

		k = &sincKernelConfig[1][i];
		makeSincKernel(fSinc16[i], 16, SINC_PHASES, k->kaiserBeta, k->sincCutoff);
	}

	// resampling ratios for sinc kernel selection
	sincRatio1 = (uint64_t)(1.1875 * MIXER_FRAC_SCALE); // fSinc[0] if <=
	sincRatio2 = (uint64_t)(1.5000 * MIXER_FRAC_SCALE); // fSinc[1] if <=, else fSinc[2] if >

	return true;
}

void freeWindowedSincTables(void)
{
	for (int32_t i = 0; i < SINC_KERNELS; i++)
	{
		if (fSinc8[i] != NULL)
		{
			free(fSinc8[i]);
			fSinc8[i] = NULL;
		}

		if (fSinc16[i] != NULL)
		{
			free(fSinc16[i]);
			fSinc16[i] = NULL;
		}
	}
}
